legend('bottom',
legend = c("Rely on informal", "Rely on formal", "Rely on UNMIL"),
bty = "o", ncol = 3, inset = 0.05, cex = 0.9, pt.cex = 1, pch = c(15, 16, 17))
dev.off()
#####################################################################################################
#
#   TITLE: International Intervention and the Rule of Law after Civil War: Replication Files (Part 2)
#   AUTHOR: Robert A. Blair
#
#####################################################################################################
# Clear
rm(list = ls())
# Set directory
#setwd("YOUR DIRECTORY")
setwd("/Users/robertblair/Dropbox (Personal)/Documents/Research/Quant/Paper/Submissions/IO Acceptance/Analysis Replication Files")
# Load packages
library(arm)
library(data.table)
library(MASS)
library(foreign)
# Replicate Figure 1
df = fread("r2_res_freqtype_pref_mlogit_margins.txt", skip = 1, header = T, check.names=T)
rows = df[c(3,5,7,9,11)]$V1
cols = c("State", "Local", "UNMIL")
r2_c_pref_all = df[c(3,5,7,9,11),
c('r2_c_pref_all','r2_c_pref_all.1','r2_c_pref_all.2')]
r2_c_pref_all = as.matrix(sapply(r2_c_pref_all, as.numeric))
rownames(r2_c_pref_all) = rows
colnames(r2_c_pref_all) = cols
r2_c_pref_all.se = df[c(4,6,8,10,12),
c('r2_c_pref_all','r2_c_pref_all.1','r2_c_pref_all.2')]
r2_c_pref_all.se = as.matrix(sapply(r2_c_pref_all.se, as.numeric))
rownames(r2_c_pref_all.se) = rows
colnames(r2_c_pref_all.se) = cols
pdf("Figure1.pdf", 10, 7)
par(family="Times")
coefplot(as.numeric(r2_c_pref_all[,2]), as.numeric(r2_c_pref_all.se[,2]), CI=2,
ylab="Change in predicted probability", main = "", mar = c(1,5.1,5.1,1), ylim=c(-.3, .3),
xlim = c(1,5.5),vertical=FALSE, pch.pts = 15, var.las = 1, cex.pts = 1.5, cex.var =1,
varnames = c('            weekly\n            patrols', '            monthly\n            patrols',
'            occasional\n            patrols', '        interventions',
'            public\n            works'))
coefplot(as.numeric(r2_c_pref_all[,1]), as.numeric(r2_c_pref_all.se[,1]),
vertical=FALSE, pch.pts = 16, cex.pts = 1.5, cex.var =1,
add=TRUE, offset = 0.1)
coefplot(as.numeric(r2_c_pref_all[,3]), as.numeric(r2_c_pref_all.se[,3]),
vertical=FALSE, pch.pts = 17, cex.pts = 1.5, cex.var =1,
add=TRUE, offset = 0.2)
legend('bottom',
legend = c("rely on informal", "rely on formal", "rely on UNMIL"),
bty = "o", ncol = 3, inset = 0.05, cex = 0.9, pt.cex = 1, pch = c(15, 16, 17))
dev.off()
#####################################################################################################
#
#   TITLE: International Intervention and the Rule of Law after Civil War: Replication Files (Part 2)
#   AUTHOR: Robert A. Blair
#
#####################################################################################################
# Clear
rm(list = ls())
# Set directory
#setwd("YOUR DIRECTORY")
setwd("/Users/robertblair/Dropbox (Personal)/Documents/Research/Quant/Paper/Submissions/IO Acceptance/Analysis Replication Files")
# Load packages
library(arm)
library(data.table)
library(MASS)
library(foreign)
# Replicate Figure 1
df = fread("r2_res_freqtype_pref_mlogit_margins.txt", skip = 1, header = T, check.names=T)
rows = df[c(3,5,7,9,11)]$V1
cols = c("State", "Local", "UNMIL")
r2_c_pref_all = df[c(3,5,7,9,11),
c('r2_c_pref_all','r2_c_pref_all.1','r2_c_pref_all.2')]
r2_c_pref_all = as.matrix(sapply(r2_c_pref_all, as.numeric))
rownames(r2_c_pref_all) = rows
colnames(r2_c_pref_all) = cols
r2_c_pref_all.se = df[c(4,6,8,10,12),
c('r2_c_pref_all','r2_c_pref_all.1','r2_c_pref_all.2')]
r2_c_pref_all.se = as.matrix(sapply(r2_c_pref_all.se, as.numeric))
rownames(r2_c_pref_all.se) = rows
colnames(r2_c_pref_all.se) = cols
pdf("Figure1.pdf", 10, 7)
par(family="Times")
coefplot(as.numeric(r2_c_pref_all[,2]), as.numeric(r2_c_pref_all.se[,2]), CI=2,
ylab="Change in predicted probability", main = "", mar = c(1,5.1,5.1,1), ylim=c(-.3, .3),
xlim = c(1,5.5),vertical=FALSE, pch.pts = 15, var.las = 1, cex.pts = 1.5, cex.var =1,
varnames = c('            weekly\n            patrols', '            monthly\n            patrols',
'            occasional\n            patrols', '        interventions',
'            public\n            works'))
coefplot(as.numeric(r2_c_pref_all[,1]), as.numeric(r2_c_pref_all.se[,1]),
vertical=FALSE, pch.pts = 16, cex.pts = 1.5, cex.var =1,
add=TRUE, offset = 0.1)
coefplot(as.numeric(r2_c_pref_all[,3]), as.numeric(r2_c_pref_all.se[,3]),
vertical=FALSE, pch.pts = 17, cex.pts = 1.5, cex.var =1,
add=TRUE, offset = 0.2)
legend('bottom',
legend = c("rely on informal", "rely on formal", "rely on UNMIL"),
bty = "o", ncol = 3, inset = 0.05, cex = 0.9, pt.cex = 1, pch = c(15, 16, 17))
dev.off()
# Replicate Figure 2
df_anyres = fread("r2_res_anyres_pref_mlogit_margins.txt", skip = 1, header = T, check.names=T)
cols = c("State", "Local", "UNMIL")
r2_c_pref_all = df_anyres[c(3),c('r2_c_pref_all','r2_c_pref_all.1','r2_c_pref_all.2')]
colnames(r2_c_pref_all) = cols
r2_c_pref_all = as.matrix(r2_c_pref_all)
r2_c_pref_all.se = df_anyres[c(4),c('r2_c_pref_all','r2_c_pref_all.1','r2_c_pref_all.2')]
colnames(r2_c_pref_all.se) = cols
r2_c_pref_all.se = as.matrix(r2_c_pref_all.se)
pdf("Figure2.pdf", 9, 5)
par(family="Times")
coefplot(as.numeric(r2_c_pref_all[,2]), as.numeric(r2_c_pref_all.se[,2]), CI=2,
ylab="Change in predicted probability", main = "", mar = c(1,5.1,5.1,1), ylim=c(-.15, .15), xlim= c(0.85,1.45),
vertical=FALSE, pch.pts = 15, var.las = 1, cex.pts = 1.5, cex.var =1,
varnames = c('', pretty.xlabels=TRUE))
coefplot(as.numeric(r2_c_pref_all[,1]), as.numeric(r2_c_pref_all.se[,1]),
vertical=FALSE, pch.pts = 16, cex.pts = 1.5, cex.var =1,
add=TRUE, offset = 0.15)
coefplot(as.numeric(r2_c_pref_all[,3]), as.numeric(r2_c_pref_all.se[,3]),
vertical=FALSE, pch.pts = 17, cex.pts = 1.5, cex.var =1,
add=TRUE, offset = 0.3)
legend('bottom',
legend = c("rely on informal", "rely on formal",   "rely on UNMIL"),
bty = "o", ncol = 3, inset = 0.005, cex = 0.9, pt.cex = 1, pch = c(15, 16, 17))
dev.off()
# Define functions to replicate Figure 3
# The following is a slightly modified version of the Corstange LISTIT function. All the changes have been marked with comments
listit <- function(Y, X, map=NULL, method="BFGS") {
## Returns a matrix of parameters with a row for each list item
## and columns for each covariate.  The parameters are placed in
## their appropriate cells, and the remaining cells are filled
## with zeros.
place.B <- function(param, map) {
k <- nrow(map)
n.x <- ncol(map)
n.param <- length(param)
B <- matrix(NA, nrow=k,
ncol=n.x)
i.param <- 1
for (i in 1:k)
for (j in 1:n.x) {
if (map[i,j] == TRUE) {
B[i,j] <- param[i.param]
i.param <- i.param + 1
}
else
B[i,j] <- 0
}
return(B)
}
neg.LL <- function(param, Y.t, Y.c, X.t, X.c, map) {
k <- nrow(map)
B <- place.B(param=param, map=map)
phi <- apply((1 + exp(-X.t %*%
t(B[2:k,])))^(-1)
, 1, sum)
p <- as.vector(((1 + exp(-X.t %*%
B[1,]))^(-1) +
phi) / k)
q <- (1 + exp(-X.c %*% t(B[2:k,])))^(-1)
list.LL <- sum(log(choose(k, Y.t)) +
Y.t * log(p) +
(k - Y.t) * log(1 - p))
off.LL <- sum(Y.c * log(q) +
(1 - Y.c) * log(1 - q))
joint.LL <- list.LL + off.LL
return(-joint.LL)
}
neg.grad <-function(param, Y.t, Y.c, X.t, X.c, map) {
k <- nrow(map)                      #number of list items
n.x <- ncol(map)                    #number of covariates
n.grad <- sum(map)                  #number of gradients to return
B <- place.B(param=param, map=map)  #rows: list item, cols: parameters
phi <- apply((1 + exp(-X.t %*%
t(B[2:k,])))^(-1)
, 1, sum)
p <- as.vector(((1 + exp(-X.t %*%
B[1,]))^(-1) +
phi) / k)
## Calculate the derivative of p with respect to the parameters
## Loop through the list items
## For each list item, get a temporary matrix of its predictors (this.X)
## Use (i.grad) and (this.param) indices to make sure the results go into
##   the appropriate columns in the (d.p) matrix, which has a number of
##   rows equal to the number of treatment observations, and a number of
##   columns equal to the number of gradients to be calculated (i.e., the
##   number of parameters we're trying to estimate).
i.grad <- 0                         #gradient index
d.p <- matrix(0, nrow=nrow(Y.t), ncol=n.grad)
for (i in 1:k) {
this.X <- X.t[,1]
this.param <- sum(map[i,])
for (j in 2:n.x)
if (map[i,j] == TRUE)
this.X <- cbind(this.X, X.t[,j]) #only the Xs for this list item
e.vec <-
as.vector(exp(-X.t %*% B[i,]) /
(k * (1 + exp(-X.t %*% B[i,]))^(2)))
d.p[,(1 + i.grad):
(i.grad + this.param)] <-
this.X * e.vec
i.grad <- i.grad + this.param
}
## Calculate the contributions to the gradient calculations from the
## off items (i.e., the non-sensitive list items).
## Loop through the off items (items 2 to k).
## For each item, get a temporary matrix of its predictors (this.X).
## Use (i.grad) and (this.param) indices to make sure the results go into
##   the appropriate column in the (off) matrix.  (Note that the (off)
##   matrix leaves a number of initial columns as zeros, corresponding to
##   the sensitive betas which do not enter these calculations.)  The
##   number of rows is the number of control group observations, and the
##   number of columns is the number of gradients to be calculated (i.e.,
##   the number of parameters we're trying to estimate).
i.grad <- sum(map[1,])              #Start after the sensitive betas
off <- matrix(0, nrow=nrow(Y.c), ncol=n.grad)
for (i in 2:k) {
this.X <- X.c[,1]
this.param <- sum(map[i,])
for (j in 2:n.x)
if (map[i,j] == TRUE)
this.X <- cbind(this.X, X.c[,j])
q <- as.vector((1 + exp(-X.c %*%
B[i,]))^(-1))
off[,(1 + i.grad):(i.grad + this.param)] <-
(Y.c[,(i - 1)] - q) * this.X
i.grad <- i.grad + this.param
}
## Calculate all the relevant components and add them together.
## (LL.t) is the component related to the treatment group.
## (LL.off) is the component related to the off items.
## (LL.joint) is the joint components.
LL.t <- as.vector(Y.t * p^(-1)) * d.p +
as.vector((k - Y.t) * (1 - p)^(-1)) * (-d.p)
LL.t <- apply(LL.t, 2, sum)
LL.off <- apply(off, 2, sum)
LL.joint <- LL.t + LL.off
return(-LL.joint)
}
## Removes observations that have missing values from the
## working data.
na.rm <- function(Y=Y, X=X) {
k <- ncol(Y)                          #number of list items
x.names <- colnames(X)                #names of the covariates
X <- cbind(1,X)                       #append a constant to X
colnames(X) <- c("constant", x.names) #names the columns
## D is a data frame.  The first column will be a vector of
## dummies indicating whether these observations will be used
## or not depending on whether or not they have NAs among the
## covariates.  The second column marks whether the observation
## is in the treatment group (ie, it gets the list) or the
## control group (it gets the off items individually).  The
## Y values follow, then the X values (with an appended constant).
D <- as.data.frame(cbind(NA,NA,Y,X))
colnames(D) <- c("to.use", "control",
colnames(Y),
colnames(X))
## If one of the off items is NA, puts NA in control, else a number
D[,"control"] <-
apply(D[,(4:(k + 2))], 1, sum)
## If control is not NA, it's a control variable, else it's NA...
D[,"control"] <-
ifelse(!is.na(D[,"control"]),
TRUE, D[,"control"])
## ...But if there is a list count, it's not a control variable
D[,"control"] <-
ifelse(!is.na(D[,3]),
FALSE, D[,"control"])
## Puts a count of the number of covariate NAs into to.use
D[,"to.use"] <-
apply(is.na(D[,((k + 3):ncol(D))]),
1, sum)
## Mark to.use if there are no covariate NAs on this obs
D[,"to.use"] <-
ifelse((D[,"to.use"] != 0),
FALSE, TRUE)
## Mark to.use as FALSE if control observations have NA Y values
D[,"to.use"] <-
ifelse(is.na(D[,"control"]),
FALSE, D[,"to.use"])
## Take the subset of usable observations, stripping off
## the to.use and control columns
D <- subset(D[3:ncol(D)], (D[,"to.use"] == TRUE))
return(D)                             #return the usable Y and X values
}
## Checks for the map defining which covariates predict which list items.
## If no map is supplied by the user, it assumes that every covariate in X
## predicts every list item.  If a map is supplied by a user, it appends
## an initial column of "TRUE" values to indicate that the constant is
## a covariate for every list item.
if (is.null(map))
map <- matrix(TRUE, nrow=ncol(Y),
ncol=(ncol(X) + 1))
else map <- cbind(TRUE, map)
## Processes the data.  First it removes observations with NA
## values, then places into separate objects control and treatment
## group observations on Y and X.
k <- ncol(Y)
YX <- na.rm(Y=Y, X=X)
YX.c <- subset(YX, is.na(YX[,1]))
YX.t <- subset(YX, !is.na(YX[,1]))
Y.c <- as.matrix(YX.c[,2:k])
X.c <- as.matrix(YX.c[(k + 1):ncol(YX.c)])
Y.t <- as.matrix(YX.t[,1])
X.t <- as.matrix(YX.t[(k + 1):ncol(YX.t)])
## The number of parameters are the sum total TRUE values
## in the parameter map
n.param <- sum(map)
result <- optim(par=rep(0, n.param),
fn=neg.LL, gr=neg.grad,
hessian=TRUE,
method=method,
Y.t=Y.t, Y.c=Y.c,
X.t=X.t, X.c=X.c,
map=map)
neg.se <- FALSE
for (i in 1:ncol(X.t))
if (is.nan(sqrt(diag(solve(result$hessian)))[i])) {
neg.se <- TRUE
break
}
output <- list(b.star=(result$par)[1:ncol(X.t)],
var=solve(result$hessian),
se=sqrt(diag(solve(result$hessian)))[1:ncol(X.t)],
vcov=solve(result$hessian)[1:ncol(X.t),1:ncol(X.t)], # added line
neg.se=neg.se,
converged=(result$convergence == 0),
n.treatment=nrow(Y.t),
n.control=nrow(Y.c),
items=k,
logL=-result$value,
deviance=2 * result$value,  ## ")" was removed and "," added to the end of this line
## the following three lines were added:
b.control=(result$par)[(ncol(X.t)+1):length(result$par)],
se.control=sqrt(diag(solve(result$hessian)))[(ncol(X.t)+1):length(result$par)],
vcov.control=solve(result$hessian)[(ncol(X.t)+1):length(result$par),(ncol(X.t)+1):length(result$par)])
class(output) <- "listit"
return(output)
}
logistic <- function(x) exp(x)/(1+exp(x))
monte.carlo<-function(freq, freq.not, n.draws, draws){
freq.list <- freq.not.list <- freq.diff <- rep(NA, n.draws)
for (d in 1:n.draws) {
par.g <- draws[d, ]
pred.freq.list <- na.exclude(logistic(freq %*% par.g))
pred.freq.not.list <- na.exclude(logistic(freq.not %*% par.g))
freq.list[d] <- mean(pred.freq.list)
freq.not.list[d] <- mean(pred.freq.not.list)
freq.diff[d] <- freq.list[d] - freq.not.list[d]
}
freq.diff
}
find.coefs<-function(sensitive, control1, control2, control3, vcov, vcov1, vcov2, vcov3, n.draws=10000){
draws <- mvrnorm(n = n.draws, mu = sensitive, Sigma = vcov) #sensitive item
freq.diff <- monte.carlo(weekly, weekly.not, n.draws, draws)
weekly.diff.mean <- mean(freq.diff)
weekly.diff.se <- sd(freq.diff)
freq.diff <- monte.carlo(monthly, monthly.not, n.draws, draws)
monthly.diff.mean <- mean(freq.diff)
monthly.diff.se <- sd(freq.diff)
freq.diff <- monte.carlo(occasional, occasional.not, n.draws, draws)
occasional.diff.mean <- mean(freq.diff)
occasional.diff.se <- sd(freq.diff)
freq.diff <- monte.carlo(interventions, interventions.not, n.draws, draws)
interventions.diff.mean <- mean(freq.diff)
interventions.diff.se <- sd(freq.diff)
freq.diff <- monte.carlo(public.works, public.works.not, n.draws, draws)
public.works.diff.mean <- mean(freq.diff)
public.works.diff.se <- sd(freq.diff)
draws <- mvrnorm(n = n.draws, mu = control1, Sigma = vcov1) #Call police
freq.diff <- monte.carlo(weekly, weekly.not, n.draws, draws)
weekly.diff.mean1 <- mean(freq.diff)
weekly.diff.se1 <- sd(freq.diff)
freq.diff <- monte.carlo(monthly, monthly.not, n.draws, draws)
monthly.diff.mean1 <- mean(freq.diff)
monthly.diff.se1 <- sd(freq.diff)
freq.diff <- monte.carlo(occasional, occasional.not, n.draws, draws)
occasional.diff.mean1 <- mean(freq.diff)
occasional.diff.se1 <- sd(freq.diff)
freq.diff <- monte.carlo(interventions, interventions.not, n.draws, draws)
interventions.diff.mean1 <- mean(freq.diff)
interventions.diff.se1 <- sd(freq.diff)
freq.diff <- monte.carlo(public.works, public.works.not, n.draws, draws)
public.works.diff.mean1 <- mean(freq.diff, na.rm=TRUE)
public.works.diff.se1 <- sd(freq.diff, na.rm=TRUE)
draws <- mvrnorm(n = n.draws, mu = control2, Sigma = vcov2) #Call UNMIL
freq.diff <- monte.carlo(weekly, weekly.not, n.draws, draws)
weekly.diff.mean2 <- mean(freq.diff)
weekly.diff.se2 <- sd(freq.diff)
freq.diff <- monte.carlo(monthly, monthly.not, n.draws, draws)
monthly.diff.mean2 <- mean(freq.diff)
monthly.diff.se2 <- sd(freq.diff)
freq.diff <- monte.carlo(occasional, occasional.not, n.draws, draws)
occasional.diff.mean2 <- mean(freq.diff)
occasional.diff.se2 <- sd(freq.diff)
freq.diff <- monte.carlo(interventions, interventions.not, n.draws, draws)
interventions.diff.mean2 <- mean(freq.diff)
interventions.diff.se2 <- sd(freq.diff)
freq.diff <- monte.carlo(public.works, public.works.not, n.draws, draws)
public.works.diff.mean2 <- mean(freq.diff)
public.works.diff.se2 <- sd(freq.diff)
draws <- mvrnorm(n = n.draws, mu = control3, Sigma = vcov3) #Do nothing
freq.diff <- monte.carlo(weekly, weekly.not, n.draws, draws)
weekly.diff.mean3 <- mean(freq.diff)
weekly.diff.se3 <- sd(freq.diff)
freq.diff <- monte.carlo(monthly, monthly.not, n.draws, draws)
monthly.diff.mean3 <- mean(freq.diff)
monthly.diff.se3 <- sd(freq.diff)
freq.diff <- monte.carlo(occasional, occasional.not, n.draws, draws)
occasional.diff.mean3 <- mean(freq.diff)
occasional.diff.se3 <- sd(freq.diff)
freq.diff <- monte.carlo(interventions, interventions.not, n.draws, draws)
interventions.diff.mean3 <- mean(freq.diff, na.rm=TRUE)
interventions.diff.se3 <- sd(freq.diff, na.rm=TRUE)
freq.diff <- monte.carlo(public.works, public.works.not, n.draws, draws)
public.works.diff.mean3 <- mean(freq.diff, na.rm=TRUE)
public.works.diff.se3 <- sd(freq.diff, na.rm=TRUE)
means <- data.frame(c(weekly.diff.mean, weekly.diff.mean1, weekly.diff.mean2, weekly.diff.mean3),
c(monthly.diff.mean, monthly.diff.mean1, monthly.diff.mean2, monthly.diff.mean3),
c(occasional.diff.mean, occasional.diff.mean1, occasional.diff.mean2, occasional.diff.mean3),
c(interventions.diff.mean, interventions.diff.mean1, interventions.diff.mean2, interventions.diff.mean3),
c(public.works.diff.mean, public.works.diff.mean1, public.works.diff.mean2, public.works.diff.mean3))
rownames(means) <- cbind('Use trial by ordeal', 'Call police', 'Call UNMIL', 'Do nothing')
colnames(means) <- c('weekly', 'monthly', 'occasional', 'interventions', 'public works')
ses <- data.frame(c(weekly.diff.se, weekly.diff.se1, weekly.diff.se2, weekly.diff.se3),
c(monthly.diff.se, monthly.diff.se1, monthly.diff.se2, monthly.diff.se3),
c(occasional.diff.se, occasional.diff.se1, occasional.diff.se2, occasional.diff.se3),
c(interventions.diff.se, interventions.diff.se1, interventions.diff.se2, interventions.diff.se3),
c(public.works.diff.se, public.works.diff.se1, public.works.diff.se2, public.works.diff.se3))
rownames(ses) <- cbind('Use trial by ordeal', 'Call police', 'Call UNMIL', 'Do nothing')
colnames(ses) <- c('weekly', 'monthly', 'occasional', 'interventions', 'public works')
list(means,ses)
}
plot.margins <- function(Y,X){
listit.model  <- listit(Y=Y, X=X)
num <- length(listit.model$b.star)
sensitive <- listit.model$b.star
control1 <- listit.model$b.control[1:num]
control2 <- listit.model$b.control[(1+num):(2*num)]
control3 <- listit.model$b.control[(1+2*num):(3*num)]
vcov <- listit.model$vcov
vcov1 <- listit.model$vcov.control[1:num,1:num]
vcov2 <- listit.model$vcov.control[(1+num):(2*num),(1+num):(2*num)]
vcov3 <- listit.model$vcov.control[(1+2*num):(3*num),(1+2*num):(3*num)]
weekly <<- cbind(1,as.matrix(subset(X, X[,1] == 1)))
weekly.not <<- cbind(1,as.matrix(subset(X, X[,1] == 0)))
monthly <<- cbind(1,as.matrix(subset(X, X[,2] == 1)))
monthly.not <<- cbind(1,as.matrix(subset(X, X[,2] == 0)))
occasional <<- cbind(1,as.matrix(subset(X, X[,3] == 1)))
occasional.not <<- cbind(1,as.matrix(subset(X, X[,3] == 0)))
interventions <<- cbind(1,as.matrix(subset(X, X[,4] == 1)))
interventions.not <<- cbind(1,as.matrix(subset(X, X[,4] == 0)))
public.works <<- cbind(1,as.matrix(subset(X, X[,5] == 1)))
public.works.not <<- cbind(1,as.matrix(subset(X, X[,5] == 0)))
vals = find.coefs(sensitive, control1, control2, control3, vcov, vcov1, vcov2, vcov3)
par(family="Times")
coefplot(t(vals[[1]])[,1], t(vals[[2]])[,1], CI=2, ylim=c(-.7, .7), xlim = c(1,5.5), vertical=FALSE, pch.pts = 15,
var.las = 1, cex.pts = 1.5, cex.var =1, mar = c(1,5.1,5.1,1), ylab="Change in predicted probability", main = " ",
varnames = c('            weekly\n            patrols', '            monthly\n            patrols',
'            occasional\n            patrols', '        interventions',
'            public\n            works'))
coefplot(t(vals[[1]])[,2], t(vals[[2]])[,2], vertical=FALSE, pch.pts = 16, cex.pts = 1.5, cex.var =1, add=TRUE, offset = 0.1)
coefplot(t(vals[[1]])[,3], t(vals[[2]])[,3], vertical=FALSE, pch.pts = 17, cex.pts = 1.5, cex.var =1, add=TRUE, offset = 0.2)
coefplot(t(vals[[1]])[,4], t(vals[[2]])[,4], vertical=FALSE, pch.pts = 18, cex.pts = 1.5, cex.var =1, add=TRUE, offset = 0.3)
legend('top',  legend = c("use trial by ordeal", "call police", "call UNMIL", "do nothing"),  bty = "o",
ncol = 4, cex = 0.7, pt.cex = 1, pch = c(15, 16, 17, 18))
}
# Replicate Figure 3
list = read.dta("survey-r3_leaders.dta", convert.factors=TRUE,convert.underscore=TRUE)
Y.busthap = cbind(list$r3.l.litrbusthap.count, list$r3.l.licobustcophap.num, list$r3.l.licobustunhap.num, list$r3.l.licobustleavehap.num)
R3.C.UNMIL = cbind(list$r3.l.cunmil.week, list$r3.l.cunmil.month, list$r3.l.cunmil.rare,list$r3.l.cunmilpala, list$r3.l.cunmilmaintain)
R2.C.CONTROLS = cbind(list$r2.l.croad.dist.rainy,   list$r2.l.cfacilities,  list$r2.l.ccapital.offense,
list$r2.l.ccoll.viol, list$r2.c.canycrime, list$r2.l.cfreq.visitpol)
colnames(R3.C.UNMIL) = c("r3.l.cunmil.week", "r3.l.cunmil.month", "r3.l.cunmil.rare3","r3.l.cunmilpala", "r3.l.cunmilmaintain")
colnames(R2.C.CONTROLS) = c("r2.l.croad.dist.rainy", "r2.l.cfacilities", "r2.l.ccapital.offense",
"r2.l.ccoll.viol", "r2.c.canycrime", "r2.l.cfreq.visitpol")
X.busthap <- cbind(R3.C.UNMIL, R2.C.CONTROLS)
pdf("Figure3.pdf", 10, 7)
plot.margins(Y.busthap, X.busthap)
dev.off()
# Replicate Figure 4
df_persist = fread("r3_res_freqtype_pref_mlogit_margins.txt", skip = 1, header = T, check.names=T)
rows = df_persist[c(3,5,7,9,11)]$V1
cols = c("State", "Local", "UNMIL")
r3_c_pref_all = df_persist[c(3,5,7,9,11),
c('r3_c_pref_all','r3_c_pref_all.1','r3_c_pref_all.2')]
r3_c_pref_all = as.matrix(sapply(r3_c_pref_all, as.numeric))
rownames(r3_c_pref_all) = rows
colnames(r3_c_pref_all) = cols
r3_c_pref_all.se = df_persist[c(4,6,8,10,12),
c('r3_c_pref_all','r3_c_pref_all.1','r3_c_pref_all.2')]
r3_c_pref_all.se = as.matrix(sapply(r3_c_pref_all.se, as.numeric))
rownames(r3_c_pref_all.se) = rows
colnames(r3_c_pref_all.se) = cols
pdf("Figure4.pdf", 10, 7)
par(family="Times")
coefplot(as.numeric(r3_c_pref_all[,2]), as.numeric(r3_c_pref_all.se[,2]), CI=2,
ylab="Change in predicted probability", main = "", mar = c(1,5.1,5.1,1), ylim=c(-.3, .3),
xlim = c(1,5.5),vertical=FALSE, pch.pts = 15, var.las = 1, cex.pts = 1.5, cex.var =1,
varnames = c('            weekly\n            patrols', '            monthly\n            patrols',
'            occasional\n            patrols', '        interventions',
'            public\n            works'))
coefplot(as.numeric(r3_c_pref_all[,1]), as.numeric(r3_c_pref_all.se[,1]),
vertical=FALSE, pch.pts = 16, cex.pts = 1.5, cex.var =1,
add=TRUE, offset = 0.1)
coefplot(as.numeric(r3_c_pref_all[,3]), as.numeric(r3_c_pref_all.se[,3]),
vertical=FALSE, pch.pts = 17, cex.pts = 1.5, cex.var =1,
add=TRUE, offset = 0.2)
legend('bottom',
legend = c("rely on informal", "rely on formal", "rely on UNMIL"),
bty = "o", ncol = 3, inset = 0.05, cex = 0.9, pt.cex = 1, pch = c(15, 16, 17))
dev.off()
